Velocity correlations in granular materials 
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A system of inelastic hard disks in a thin pipe capped by hot walls is studied with the aim of 
investigating velocity correlations between particles. Two effects lead to such correlations: inelastic 
collisions help to build localized correlations, while momentum conservation and diffusion produce 
long ranged correlations. In the quasi-elastic limit, the velocity correlation is weak, but it is still 
important since it is of the same order as the deviation from uniformity. For system with stronger 
inelasticity, the pipe contains a clump of particles in highly correlated motion. A theory with 
empirical parameters is developed. This theory is composed of equations similar to the usual hydro- 
dynamic laws of conservation of particles, energy, and momentum. Numerical results show that the 
theory describes the dynamics satisfactorily in the quasi-elastic limit, however only qualitatively for 
stronger inelasticity. 



PACS numbers; 81.05.Rm, 05.20.Dd, 47.50. +d 



I. INTRODUCTION 



A granular system normally consists of a large number 
of particles colliding with one another and losing a lit- 
tle energy in each collision. If such a system is shaken to 
keep it in motion, its dynamics resembles that of fluids, in 
that the grains move seemingly randomly. Considerable 
effort has been devoted to the development of a contin- 
uum description for hydrodynamic equations [p]-pl|. 

Two approaches are employed by different authors. 
One is to set up a Boltzmann equation and then 

to calculate hydrodynamic quantities by doing averag- 
ing with the distribution function derived from the equa- 
tion. In this case, the molecular chaos assumption of 
the Boltzmann equation assumes zero correlations be- 
tween particles. The other approach is to specify some 
hydrodynamic quantities, and then write down the con- 
servation equations for them [p|-|rT| . Generally, there are 
three equations: conservation of mass, balance of energy, 
conservation of momentum. The mass conservation is in 
the standard hydrodynamic form. The momentum flux 
balance equation is in the form of Navier-Stokes equation 
for fluid dynamics. The energy "conservation" equation 
includes dissipation of energy via collisions. 

The failing of Liouville's theorem for granular systems 
I p^ casts doubts on the applicability of conventional ap- 
proaches to hydrodynamic equations. In stead of writing 
down such equations based on unjustified assumptions, 
studying the dynamics with as few assumptions as pos- 
sible, and trying to develop a theory closely connected 
to experiments, may be a less ambitious, but more solid 
approach. 

One major consequence of the usual hydrodynamic 
theories of fluids is the Maxwell-Boltzmann distribution. 
In the frame in which the average system velocity is 
zero, this distribution implies no momentum correlations 
whatsoever among different particles. This result is true 
for classical particles independent of the strength of the 
inter-particle potential. In contrast, however, granular 



materials commonly induce correlated collective behav- 
iors. Think about the surface waves of vibrated sand , 
or their convection patterns [ p^ , (for a recent review, see 
p^.) The grains which take part in these collective be- 
haviors do have correlated velocities. We therefore ask: 
how important are these correlations and how are they 
built up? 

In this paper, we investigate the building up of correla- 
tions between velocities of grains. There are two mecha- 
nisms upon which correlations can be built up. One is the 
inelastic collisions between particles — after a collision, 
the velocity difference between two particles is smaller 
than that before the collision. This is a local effect, and 
the correlation is short ranged. The other mechanism is 
from momentum conservation — the larger the scale of a 
perturbation producing a localized average velocity, the 
slower the perturbation decays p^ . Fluctuations make 
the system non-uniform, so that there are localized clus- 
ters of particles all moving with about the same velocity. 
This correlation effect is a result of fluctuations which 
are neglected in the usual hydrodynamic treatments. 

There are hydrodynamic theories ignoring fluctuations 
which are consistent with numerical results for weak in- 
elasticities, but are quite inaccurate when inelasticities 
are strong, see, e.g. [0. These theories are attempts 
to describe velocity fluctuations about some mean flow. 
They work flne in quasi-elastic regime because correla- 
tions are small and negligible. But when inelasticities are 
strong, collisions can bring groups of neighboring parti- 
cles to essentially the same velocity and thereby produce 
a correlated motion which enhances the observable ef- 
fect of any fluctuations in the system. We shall see this 
happen in our study. 

The boundary conditions and system sizes indepen- 
dence of the essential characteristics of thermodynamic 
systems is one indication that thermodynamics is a uni- 
versal description. However, this independence is lost in 
granular systems. We show that a universal description 
may still exist for the unconserved modes of the dynam- 
ics. 
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II. THE THIN PIPE MODEL 

A. The system 

To investigate the validity of a hydrodynamic descrip- 
tion, we should study the simplest situation which can 
show hydrodynamic behavior. In the elastic case, one di- 
mensional systems have too many conservation laws and 
do not show a fully ergodic or hydrodynamic behavior 
[^,|8|. Here, we investigate a two-dimensional system in 
the form of a long thin pipe (Fig. T]). The grains confined 
in the pipe are all identical, and the width of the pipe 
is set so that two grains cannot pass each other. Thus 
the motion of grains is two dimensional to ensure ergod- 
icity, while at the same time we can order these grains. 
Pipe systems were studied before However, the no- 
passing condition enable this thin pipe model to simplify 
greatly both numerical and analytical calculations. (A 
two-dimensional version of this model is studied by Brey 
and Cubero |].) 



FIG. 1. Snapshots of the thin pipe system. The periodic 
side walls are indicated by dashed lines. The two end walls 
are energy sources kept at the same temperature. Notice how 
most of the particles fall into a cluster, which moves up and 
back through the system. 

The two side walls are periodic — after leaving one side 
wall a particle comes back through the other. The dis- 
tance between the side walls is chosen to be 2.5 times the 
radius of a particle. This choice prevents any passing. 
Two end walls are energy sources, and are kept at the 
same temperature. 

For a thermodynamic system, the bulk properties 
should not depend on the details of boundary conditions. 
However, for some granular properties, boundary con- 
ditions can be quite important We employ two 
different boundary conditions in the numerical calcula- 
tions: In both cases, when a particle hits an end wall the 



direction of its motion is turned around, and the parti- 
cle is returned to the system. In the fixed speed bound- 
ary condition, the returned particles move away from the 
wall with a unit speed. Alternatively, in the Boltzmann 
boundary condition the returning speed is picked from a 
distribution P{u) = 2ue-"' |]. All fi gures describe sim- 
ulations with the Boltzmann boundary conditions unless 
otherwise specified. 



B. The parameters and variables 

We use the simplest model: non-rotating particles. Af- 
ter a collision, the radial relative velocity changes sign, 
and decreases by a factor of the restitution coefficient r, 
with < r < 1. In the collision, the other components 
of the velocities are unchanged. Thus, r = 1 is for elastic 
particles, and r = for extremely inelastic particles. We 
also define e = 1 — r. 

The coordinate in the problem is an index, i, which 
indicates the position of the particle. Suppose there are 
N = 2n particles in the thin pipe. They are ordered as 



-n,-{n- 1), 



-l,0,I,---,(n-l). 



By using the particle number, i, as our coordinate, we 
take advantage of the 'no-passing' property of the thin 
pipe and thereby get a Lagrangian description of the sys- 
tem. 

Let us denote the velocity of the ith particle as Ui, the 
relative velocity between the ith and the (i-l-l)th particle 
as Vi = Ui+i — Ui, and the velocity of the center of mass 
as u, the velocity of the ith particle with respect to the 
center of mass as it^ = Ui — u. Let us assume the pipe 
is along the x direction. Then the x-component of the 
velocities are special and we denote them as Ui and Vi. 
We use an over-line notation to indicate the root mean 
square (rms) value of some quantities. 

We propose a method to calculate profiles of various 
quantities throughout the system and the velocity corre- 
lations. (A profile is a plot of the value of some averaged 
quantity as a function of the particle number variable 
i.) Instead of the strongly correlated velocities m's, we 
study the relative velocities of neighboring particles, i/. 
In using w we focus our attention on the relative motion 
of the particles and away from their collective and corre- 
lated motion. 

There are four parameters which will describe our sys- 
tem, the particle number, N, the pipe length L, width 
W, and the inelasticity e. Of course e measures the total 
amount of inelasticity in one collision. In a system with 
many particles, the effect of the inelasticity is enhanced 
by the correlation effects. For this reason, we expect two 
combinations of N and e to be important. The product 
Ne measures the total amount of inelasticity in the sys- 
tem. For a one dimensional system, imagine a particle 
with a large velocity hitting a group of n particles, sitting 
almost at rest. The added momentum will cascade down 
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the group until at the end of the Hne the transmitted 
momentum will be diminished by a factor exp(— ne). In 
addition, a previous calculation showed that dissipa- 
tion of energy led to a gradual decay of temperature in 
the form of an exponential of —c^n where c is a con- 
stant. Thus we expect a dip in temperature determined 
by the combination of parameters, \flN . Changing the 
remaining parameters, L and W ^ will only modify some 
numerical factors in the theory — but will not change the 
qualitative behavior of the system. 

The system showed in Figure |l| contains both low den- 
sity and intermediate density regimes. There are some 
complication in such systems because of different geo- 
metrical factors for different density regimes To 
avoid such complication and focus on the dynamics of 
the system, we carry out our numerical calculations only 
for systems with extremely high density, where the typ- 
ical spacing between neighboring particles is about 2% 
of the radius of a particle; or for systems with extremely 
low density, where the spacing at the highest density re- 
gion of the system is about 10 times the radius. The 
essential characteristics of the dynamics are independent 
of density regimes. 



C. The steady state 

This system can reach a statistical steady state. In this 
state, the particles move fast near the hot walls, and the 
density is low there. Towards the center of the system, 
the density is higher. For quasi-elastic situations, the sys- 
tem is relatively uniform; but for stronger inelasticities, 
the particles near the center can form a cluster and move 
with about the same velocity. The cluster was seen and 
understood in previous calculations | pT| , p^ . The relative 
motion of particles is reduced by the inelastic collisions 
between them. In fact, when Ne is large, the relative 
motion can be very small and then momentum conser- 
vation causes that each particle in the cluster has about 
the same velocity which is just the mean velocity of the 
cluster. 
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FIG. 2. Profiles of a low density system of 100 particles 
with r — 0.94, just above the critical value for inelastic col- 
lapse rc. o is for (tti), * for {v1/2), and + for {mui+i). 

Figure ^ shows a plot of some profiles in an inelastic 
situation with two hot walls. Notice that the profile of 
(wf ) has a flat region at the center. This was seen before 
Q. That flattening occurs because the central particles 
almost always fall within a cluster, and the cluster moves 
with a large average velocity but small relative velocities. 
The plot of (w^) indeed shows that the relative velocity 
decreases to a very small value near the center of the sys- 
tem. This decay in (vf) is roughly what we might expect 
from a simple hydrodynamic description, in which one 
balances energy flux with dissipation [0 . The hydrody- 
namics then gives an i-dependence which is a superposi- 
tion of growing and decaying exponentials. That theory 
is in some sense a mean field theory which ignores the 
correlations between velocities. In much of what we do, 
delicate and long-range correlations effects will be very 
important for the behavior of u's but less important for 
the v's. In fact we shall see that the rms of vt obeys 

—Vi = b'^v,, (1) 



where is proportional to e for small values of the in- 
elasticity. The solution to the equation is 

Vi = vq cosh(6i). (2) 

Equations (|l]) and (Q) describe a situation in which heat 
conduction balances against energy dissipation. 

On the other hand, the large degree of correlation be- 
tween Ui and Mi+i is quite unexpected. No such correla- 
tion occurs in the usual statistical mechanics. This kind 
of correlation effect is not directly contained in any hy- 
drodynamic equations. As we shall see, it is a result of 
fluctuations not usually included in hydrodynamics. 
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FIG. 3. Profiles of {v1/2) and {u}/2) for two diflterent 
boundary conditions. The system has low density 100 par- 
ticles and r = 0.94. Each profile is rescaled by changing the 
scale of velocity so that {vq/2) = 1. There are two lines, 
which nearly overlap each other, describing (ii|/2). The o is 
for (m?) with the Boltzmann boundary condition, and the + 
is for {u1) with fixed speed boundary condition. These two 
profiles are very different. 



indicates the amount of correlation in the motion. When 
the inelasticity is weak, the velocity correlations are also 
weak, and this ratio is very close to unity. For strong 
inelasticity, where correlations are strong, this ratio can 
be very small. 



Boundary conditions are often important for granular 
systems. Figure || shows the effects of boundary con- 
ditions. We see (uf) depends sensitively on boundary 
conditions. In contrast, after a rescaling, {vf) is nearly 

independent of boundary conditions. There is no similar j3_ PDF's of velociti 

rescaling which can make the profiles for {uf) overlap. 



D. Correlated motion and random motion 



Since 



= + - 2(ui-u,+i), (3) 

when the correlation (uiUi^i) is weak, we simply have 
= 2uf, assuming a weak i-dependence. But when cor- 
relation is strong, the relation between and uf is quite 
different. We shall study that difference in detail. From 
the mechanism described above, we know near the center, 
uf is roughly constant, independent oii, as a consequence 
of the motion of the cluster. Conversely, vf will vary be- 
cause of energy dissipation. In our considerations, we 
shall focus upon vf, which has an average which can be 
interpreted as a local temperature. We argue that vf is a 
more relevant variable than uf, since to a large extent, it 
determines the collision rate, and the effect of a collision. 
In addition, Vi behaves as predicted by the simple hy- 
drodynamics theory, it decays exponentially, and forms 
a hyperbolic cosine curve as a function of i. Conversely, 
Ui is produced by subtle correlation effects. 
We can also write (||) in the form 

I {uf + w'+i> = ^ {vf) + {u^u^+i). (4) 

The term on the left hand side describes the total mo- 
tion with respect to the lab frame, the second term on 
the right hand side describes the correlated motion be- 
tween particle i and particle i -I- 1, and the first term on 
the right hand side describe the random relative motion 
between neighboring particles. So put this into words, 

/ total \ _ / random ^ _|_ correlated \ 
\ motion J \ motion J \ motion J 

The first term on the right can be interpreted as a tem- 
perature; the second as a result of the correlated motion 
of the two particles. In this way, we see that the ratio 

^ {vf) random motion 

* {uf + total motion 



The probability distribution functions (PDF) for 
and Vi provide considerable additional insight into the 
nature of the system. See figures (Fig. ^) and (Fig. ||). 
In these figures, the variables are normalized to give each 
PDF the same variance. 

In the PDF plots for iti, we see a fundamentally Gaus- 
sian behavior inside the cluster. Outside the cluster, the 
part of the curve shown is Gaussian but there is a strong 
high velocity tail. 
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FIG. 4. Time-averages probability for Ui in a low density 
system with A'' = 100, r = 0.94. Five such curves are shown, 
which are for i = 0, —10, —20, —35, —45. The first three in 
this list are close to Gaussian and they lie almost on top of 
one another. The other two, are quite different. 



In contrast, the PDF plots for Vi show a structure 
which is essentially the same inside and outside the clus- 
ter. Thus, all over the system, the v' s behave in the same 
way, but this behavior is quite non-trivial. 
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FIG. 5. Probability distribution for relative velocities. The 
calculation is done as a time average for a low density system 
withAf = 100, r = 0.94 and i = 0, -10, -20, -35, -45. the 
PDF's for Vi's collapse into a single curve after a rescaling. 
Once again the three curves for particles inside the cluster fall 
on top of one another, while the others are slightly different. 

We will use the constancy of the PDF of Vi (in the 
whole system) and Ui (in the interior of the system) to 
develop our theoretical model. 



III. MOTION IN THE CENTER OF MASS 
FRAME 

Figure || suggests that we can decompose the dynam- 
ics of the system into two parts: I) the motion of grains 
in the center of mass frame and II) the motion of the 
center of mass itself. Part I is independent of boundary 
conditions and all the effects of boundary conditions are 
attributed to Part II. Part I is described in terms of the 
variables Vi which may be considered to be weakly corre- 
lated with one another. Part II involves variables Ui, and 
strong correlations among the variables. In this section, 
we focus our attention upon the effects of conservation 
laws upon the system, and particularly on motion of Part 
I. 



A. Theoretical calculation 

Since the number of degrees of freedom of Part I is 
equal to the number of tJ^'s, this part of motion can be 
described in terms of Vi's. So the problem can be solved 
in two steps: the rms of w^'s and the correlations between 
z/i's. Our interests in the variable t/^'s are also based on 
the numerical results showed in the previous section that 
the profile of Vi is, in accordance to hydrodynamics the- 
ory and Equation (^, a hyperbolic cosine function of i, 
plus its weak dependence on the boundary conditions — 
these suggest that Vi can form the basis of a solution to 
some simple hydrodynamics equations. 



1. Profile of Vi 

Collisions 

For the steady state, mass conservation is reduced to 
a trivial statement that (ui) = and (vi) = 0. The mo- 
mentum and energy transfer between particles are results 
of collisions between them. So to investigate momentum 
and energy conservation, we study the effects of a single 
collision first. 

Let us consider a collision between the ith and the 
{i + l)th particle, during which Ui, Ui+i, and Vi change 
to u'^, and v'^ respectively. According to the inelastic 
collision rule, 

Ui^Ui^ —nv^^n, (6) 

1 + r 

u'i+l = Ui+i —hv.^ri, (7) 

ifi = Vi - {1 + r)nvi^ri, (8) 

where n denotes a unit vector, pointing in the direction 
of the line of centers at the point of collision while Vi^n is 
the component of Vi in that direction. 
Pressure 

The collision described above results in a change in the 
momentum of particle i -\- I, 

^, ^ 1 + J- . 
Uj+i - -Uj+i = —nvi^n- 

In a long time interval t, the momentum change of par- 
ticle i + 1 from collisions between particle i and particle 
« -I- 1 is 

1 + r 

PiWt^ —Y^{n ■ x)v,^r., (9) 

where ' ') i^ summation over all the collisions be- 

tween the ith particle and the (i -I- l)th particle. The 
X component in Equation (^ is the direction along the 
pipe. 

In writing Equation (H) we have identified the rate of 
momentum transfer from particle i to particle z + 1 as 
an average pressure. Pi times the pipe-width, W ^ while 
t is the time for the summation. We shall be dealing a 
lot with sums over collisions as in equation (^). To un- 
derstand them, we should realize that ' ') / can be 
written as the rate of collisions between i and i -I- 1, q, 
times an average over collisions (• • ■)i of this type. No- 
tice that the average over collisions is very different from 
the time-average (•■■). For example, {vi) must be zero 
in any steady state situation. However, since Vi must 
be negative for a collision to occur, then {vi)i must be 
negative. 

Now go back to Equation (^). For the steady state, 
the momentum flux must be a constant, so this summa- 
tion over a long time interval must be independent of 
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i. Thus the momentum conservation law has the conse- 
quence that the pressure as defined by Equation (^) is 
independent of i. So the equation for momentum conser- 
vation in our system takes the form 



1 



r \ ^ 



2Wt 



(n • x)Vi^n 



1 



2W 



■Ci{{n ■ x)vi^n) 



p. 



(10) 



The distribution functions for relative velocity only de- 
pend weakly on i (Fig. Thus, all components and 
averages of Vi vary in proportion to one another as i is 
varied. As i approaches the center of the system, the 
typical value of the momentum transfer per collision de- 
clines in proportion to Vi. Then, by Equation (|l^) the 
collision rate increases by going inversely as the relative 
velocity. 

Using the same arguments, we can also understand 
the pressure-definition. Equation ([To|), in a familiar form. 
Pressure is a flow of momentum per unity area per unit 
time. One kind of flow involves transfer of momentum 
from the zth particle to the next one. The collision rate is 
of the order oivi/li^ where U is the mean spacing between 
the two particles. During each collision, the average mo- 
mentum transfer is proportional to TJ.;. From these two 
facts, the momentum flux is proportional to vf/h. The 
average of relative velocity squared is the temperature, 
T, while l/{Wli) is the density, p. This result is then 
in the familiar form, P — pT . This identification is an 
order of magnitude argument. For calculations, we use 
the exact result. Equation (p^. 

Energy Balance 

Now let us study the effects of this collision on the 
energy balance. The energy transfer to particle i and 
particle i -f 1 are, respectively. 



_ ^2 _ 



1 + r 
2 

1 + r 



n) 



' 2 2 



i+l,n> 



and the energy dissipation is 



We can form an energy conservation equation by bal- 
ancing the energy dissipation with the difference of the 
energy flux. However, the above expressions involve u^'s 
which are correlated and do not belong to the motion in 
the center of mass frame. To find a consistent descrip- 
tion, we want to express this conservation in terms of ViS. 
Because the essential dynamic process is determined by 
the collision rule. Equations (^[]|), an equation describing 
the balance of a quadratic form of Vi& will incorporate 
the energy conservation. 

In fact, we have, from (|[S), 



->I2 ->1 /I 2\ 2 

-^i = -(1-7- )Vi^n, 

^+1 - ^i+l = (1 + r)'Vi,nVi+l.n + 
l^-l - = (1 + r)Vi^nVt-l,n + 

4: 

For the steady state, the total change in should vanish. 



(1 + r)^, 



(1 



(i+i) a-i) 



or equivalently, 



(i+l) 



El , 1 + 2 \ 



(i-1) 



4 

1 + r 



(11) 



The first term is from energy dissipation, while the other 
terms take the form of energy transfer. 
Profile of Vi 

We wish to simplify our energy-flow equation by reduc- 
ing it to an equation for {vl^i. However, correlations 
between vi and Vi^\ appears in (p+j). We must eliminate 
these terms. For an elastic uniform system, this correla- 
tion takes a simple form. 



(UiUi+l) = ((Ui+l - Ui)(Uj+2 - Uj+l)) 



^1+1/ 



(12) 



In the elastic case, it is equally true for the usual time- 
weighted average or for the collision weighted average, 
as 



^/C CC J c c c \ 



(13) 



where vf = y X)^*' ^n/''^i^ ^^'^ total number of 

collisions between particle i and particle i + 1, = at. 
As defined here, is an collision average of Vi just before 
collisions. 

Equation ( p^ has scalars on the left and right hand 
side. There are corrections to this relation for inelastic 
particles and when there is a spatial variation in the av- 
erages. The corrections must be scalars and of order Vi'^. 
One correction is of order of the order of evi'^. In the 
other correction, ^ is applied to Vi^ . However, in virtue 
of the result in Equation (|^) these two terms are really 
the same. Consequently, we need only one of these two 
corrections. We write the resulting structure in an even 
parity form as 
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(14) 



Now Equation (|lj) can be used to transform equation 
( pi] ) into the form 



1 — flie 



or, since n^w^ is a constant independent of i, (see Equa- 
tion (p^)), we find a heat flow equation 

(6 - 4:ai)-^—Vi = Vi+i - 2vi + Vi^i. 
1 + r 

In writing the last structure we have noticed that the 
different kinds of coUision averages all have the same i- 
dependence. Now we can phrase our result in a contin- 
uum form 

(6-4ar)^.. = ^. 

In this way, we obtain 

Vi =vocosh.{bi). (15) 

where 

&2 = (3-2ai)e. (16) 

2. Correlations between velocities 

Correlations between Vi^s are short ranged. Let us only 
consider the nearest neighbor correlation. When there is 
no dissipation, the only non-vanishing correlation of the 
Vi's is the nearest neighbor average of Equation (p^. For 
r < 1, there is a small correction to that relation. Just as 
before, (see Equation (p^), we write an equation for the 
average of a nearest neighbor product in the same form 
as in the elastic case, but with a correction proportional 
to e: 



1 - 026. 



JiVi+l, 



(17) 



where the averages are time average. 

This assumption, with the profile of Vi determined 
above, completes a description of the motion of grains 
in the center of mass of frame, i.e. Part I of the dynam- 
ics described before. As an example, let us calculate the 
correlations between m"s, the velocities of particles in 
the center of mass frame. To illustrate the similarity be- 
tween this part of the dynamics and conventional thermo- 
dynamics, i.e., the independence of boundary conditions 



and system sizes, we consider the center of mass frame 
of the 2m particles at the center of the system. Keep 
in mind that rather than fixed, m can be treated as a 
variable in the following calculation. 
Express in terms of Vi's, 



2m 



^(m-j>j- {m + j)vj 

j = -{m-l) 



So 



3=1 



2«+i ^ ^ 



where 



A: 



rn — 1 



^ (m - j)vj ~ ^ (m -f j)vj 

j=i+i j = -(m — 1) 



Let us calculate the correlation between Uq and . Only 
keeping the correlations between nearest neighbor, we 
have 

m-l -1 

^ (m - j)vj - ^ (m j)vj 

_j=l j=-(m-l) 

(m — 1 m — 2 \ 

^ (m - jfv] + 2 ^ (m - j){m - j - l)v,Vj+i 

m-2 



+ Y[('m~ j)vj - (m -j - l)vj+i]^' 



+2a2£ ^ (m - j){m - j - l)vjVj+i 



So 



{2ulul) 



Ai — m? 
Ai + m? 



(18) 



Obviously, when m = 1, the above ratio is —1, because 
Mq = —u\ at all time. For elastic particles, the ratio can 
be calculated analytically to be — 2m-i ■ The inelastic- 
ity changes this dependence. Let us call 2m the 'cluster 
size', since it corresponds to the usual practice of defin- 
ing a 'cluster' then separate the motion of particles into 
mean flow and fluctuations. 

From expression (|8|), we see that when the cluster is 
large enough, Ai can be big comparing to m^, then the 
correlations between velocity fluctuations can be big. 
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B. Numerical results 



We carry out numerical simulations to investigate the 
statistical steady state of the system. Here we compare 
the numerical results with the above theory describing 
the motion in the center of mass frame. 



1. Quasi-elastic situations 



First let us look at the quasi-elastic situations, i.e. very 
small e. Before testing the profile of w^'s, we exam the 
crucial assumption, (14). 
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FIG. 6. Numerical results for 1.5 — ai for a high density 
system with 100 particles averaged over 2 x 10^ collisions, oi 
is defined in Equation (pl|). The o is for r — 0.995; the + is 
for r = 0.95. The line is for 1.5 - ai = 0.029 from the fit in 
Figure ^ 

Now let us look more sharply at the data. To find 
oi, we do a very accurate determination of the ratio of 
averages from the left and the right hand sides of equa- 
tion (p^). This equation is then solved at each z-value to 
find a local value of oi. The result is shown in Figure ||. 
The theory is right if ai is independent of i and wrong 
if it has an important i-dependence. The figure seems to 
show that there is an excellent fit for the smaller value of 
e, and a bad fit for the larger. 

From equation ( ]l6| ) we see that the important combi- 
nation determining the properties of the profile of v is 
3 — 2ai. But ai is very close to 1.5, as shown in FigurcT 
Then the ai effect changes the prefactor in equation (|lf 
from 3 to 3 — 2ai, i.e. by a factor of 50. The velocity 
correlations renormalize e, and reduce the energy dissi- 
pation. 

Also ai is essentially a local correlation effect origi- 
nated from the inelastic collisions. For an elastic system 
with comparable inhomogeneity, there is also a correction 
to the factor — ^ in Equation (|lj), but the correction is 
usually an order of magnitude smaller than the effects we 
are seeing here. 
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FIG. 7. Fit to a hyperbolic cosine curve of the profile of the 
Vi for a high density system of 100 particles and r = 0.9995. 
To check equation (|lq), we take the inverse of the hyperbolic 
cosine of Vi/vo and plot the result as a function of i. The 
straight line indicates a fit to the theory. In the theory, the 
slope is proportional to the square root of e. Here the slope 
is 0.0055, which is equal to the square root of 0.06e. 

A test of ( p^ is showed in Figure |^. Analysis like this 
permits the determination of the slope like the one in 
Figure 7 as a function of e. We have called this slope b. 
Figure S shows that the numerical values give an e de- 
pendence for b which fully supports the theory. However, 
notice that all this analysis applies to very small values 
of e. The next section considers more inelastic situations. 
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FIG. 8. The e dependence of b for systems with A'^ = 100. 
The curve is the theoretical fit, the square root of 0.058e. (See 
Equation (|l|).) 



2. Stronger inelasticity regime 

We look at smaller r's. To avoid inelastic collapses, we 
limit our r to be greater than Vc- For a system of 100 
particles with extremely high density, rc ~ 0.95. 

When r gets smaller, there is a cluster of particles mov- 
ing around the center of the pipe, all with about the 
same velocity. The system is in a state far away from 
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equilibrium. Also, it is very nonuniform — the particles 
around the center are highly correlated while those near 
the boundaries move independently; the energy flux is 
strong near the end walls, but rather weak inside the sys- 
tem. As a consequence, the PDF's of quantities change 
significantly from particles near the center to those near 
the boundaries, e.g. the PDF's of Ui's, though there is no 
big change in the PDF's of u^'s. 




particle I 

FIG. 9. Fit to a hyperbolic cosine curve of the profile of Vi 
for a high density system of 100 particles and r = 0.95. This 
is a higher-e analog of Figure ^. The straight line corresponds 
to a hyperbolic cosine profile-curve, and its slope is 0.054, a 
value extrapolated from the expression for quasi-elastic cases 
(Fig. Ifowever, the straight-line fit is not very good, espe- 
cially near the boundary. 

Figure ^ once again plots a quantity which should be 
linear in i if the theory, equations (|l^) , is right. Now, 
for this larger values of e, there are substantial varia- 
tions in slope. It appears that the theory does not apply 
for the fifteen particles nearest to each of the boundaries 
and that it might have small troubles elsewhere. This 
discrepancy is also shown when we plot the slope, cal- 
culated from doing numerical derivatives on Figure plto 
give 6 as a function of i . This plot is given as Figure^lC|. 
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FIG. 10. The position dependence of the prefactor 6 in a 
high density system with 100 particles and r — 0.95. The line 
is fe = 0.054 extrapolated from the quasi-elastic cases. 

The discrepancy between the theory and numerical re- 
sults for strong inelasticity is not surprising. Though 
taking into account the correlations between fluctua- 
tions, the theory is still based on concepts of conventional 
fluids — no internal structures are considered. However, 
when inelasticity is strong, the dynamics is affected by 
intrinsic structures of the collection of the particles, and 
the whole system may belong a different phase A 
satisfactory theory must incorporate this feature. 

Now let us look at the velocity correlations. Only the 
nearest neighbor correlation (Equation (O)) is consid- 
ered. The theory leads to the expression p8| ) of the cor- 
relation between Uq and Ui, which is independent of sys- 
tem sizes or boundary conditions. To test this expression, 
we calculate numerically this correlation with respect to 
different cluster sizes, i.e. different m, with ( [T^ ) and with 
the profile of Vi calculated numerically. The comparison 
between theory and numerical result is showed in Fig- 
ure 1 1 , We see the correlation increases with increasing 
cluster size. The comparison is the best for 02 = 0.6. 
When the cluster size is big enough, most part of the 
total motion belongs to the correlated motion. We want 
to point out that this curve is independent of bound- 
ary conditions. Also for systems with different sizes, we 
get sections of different length from this same curve, as 
shown in the figure. 
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FIG. 11. The cluster size dependence of the ratio (|l^) . The 
system is in a low density regime, with r = 0.94. * is from 
time average results of a simulation with 100 particles and o 
is from a simulation with 60 particles, and the curve is from 
([|) with aa = 0.6. 

We want to point out that the major point of Figure pd] 
is to demonstrate that part of the dynamics, the motion 
in the center of mass frame, is independent of boundary 
conditions and system sizes. The agreement between the- 
ory and numerical results can not be viewed as a strong 
support for the details of the theory because the profile of 
Vi is from numerical calculations, rather than (|l5|), also 
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the value 02 = 0.6 is a fitting parameter. The theory 
captures some qualitative features of the dynamics, but 
is still incomplete. 



IV. MOTION OF THE CENTER OF MASS 

Because the total momentum of the system can be only 
changed by the collisions between the outermost particles 
and the walls, and the motion of the outermost particles 
is close to that of a elastic system, the motion of the 
center of mass should also be close to that of a elastic 
system. For a elastic system, 



(19) 



where u* is the rms speed of the outermost particle. 
From Figure ^ we see this estimate is about right, 
though the numerical factor must be calculated from de- 
tailed distributions. The result also seems sensitive to e. 
This is because the PDF for the velocity of the outermost 
particle is more skewed for higher value of e, and so the 
ratio between u* and the momentum transfcrcd into the 
system from the wall depends on e. 

Notice that the motion of the center of mass depends 
strongly on the boundary condition. 
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FIG. 12. Test of (^9|) for two boundary conditions for high 
density systems with A'^ = 100. The ratios are all around 1, 
as we expect from our order of magnitude argument. The 
* is for the Boltzmann boundary condition, and the o is for 
the fixed speed condition. The motion of the center of mass 
depends strongly on the boundary conditions. 

Suppose the motion of particles in the center of mass 
frame is independent of the motion of the center of mass 
itself, i.e. u is uncorrelated to w^'s, then 



Simulations show that the profile of Vi is nearly indepen- 
dent of boundary conditions, so is the motion of the sys- 
tem in the center of mass frame. However, (m^) depends 
sensitively on the boundary conditions, and so does the 



motion of the particles in the lab frame, i.e. the profile 
of {uD (Fig. I). 

Due to the motion of the center of mass, the corre- 
lations between uts are enhanced, comparing to those 
between m^s. 



{2uiU. 



i+ll 



2(«+i) + 2(^2) 



(20) 



The ratio, Ri, between random motion and total motion 
was defined by us in Equation (^. 



Behavior of the ratio Ro 



When ne is small, we can expand expression (|2(]|), us- 
ing Equations (|l5|), (|l8|) and (|l^). Keeping terms linear 
in e, we have, 

- log(i?o) = ^{ [{n - + {n~ 1)] (3 - 2ai) 

+2(n- l)(n-2)a2/3}, (21) 

where < a < 1. From Equation ( pT]) we see that when 
ne is small, — log(i?o) is proportional to ne. 

Numerical results of — log(i?o) are showed in Figure [l^. 
We do see that — log(i?o) is proportional to e for very 
small e. However, when e is big, where we expect strong 
nonlinear effects, it is proportional to e^. 
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FIG. 13. The logarithm of the ratio (||) for i = ver- 
sus e. -I- is for A'^ = 100, * is for N = 70, and o is for 
A'' = 40. AH three are for low density systems with the Boltz- 
mann boundary condition. The dashed line indicates a depen- 
dence log(i?o) cx: e and the dotted line indicates a dependence 
log(i?o) oc e^. 
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As we argued in Section II, there are two important 
combinations of N and e. The product N^/e describes 
how temperature decays towards the center of the sys- 
tem, which agrees exceUently with the numerical resuhs 
when e is very small. However, Equation ( |2l| ) shows that 
in this limit, only the product Ne appears in the final 
expression for Rq. This seems to suggest that Rq, i.e., 
the degree of the coherence of the particles' motion, is 
determined by the product Ne (Fig. ^ and Fig. 15). 
These two figures exhibit rather interesting features of 
the dynamics Q , though we do not have a satisfactory 
understanding of them. 
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FIG. 14. The curves shown in Figure |l^ can be shifted to 
overlap by changing the a;-axis from e to ef{N), where f{N) 
is a function of the total number of particles in the system. 
/(lOO) = 1. Three curves are shown, they are all for low den- 
sity systems with Boltzmann boundary conditions. -I- is for 
iV = 100, * is for iV = 70 and o is for iV = 40. 
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FIG. 15. The function f{N) from Figure 0. f{N) is pro- 
portional to A'^ with a small adjustment due to boundary ef- 
fects. 



tween velocities of granular particles are shown to be im- 
portant for a correct understanding of such systems. For 
systems in quasi-elastic regime, correlation is small, but 
not negligible because the deviation from equilibrium is 
also small. For systems with stronger inelasticity, corre- 
lation is crucial for a correct theory. Our theory describes 
the dynamics satisfactorily in the quasi-elastic limit. For 
stronger inelasticities, numerical results show quite inter- 
esting behaviors of the system, however, our theoretical 
understanding is only qualitative at this stage. 

Characteristicly for granular systems, fluctuations are 
important at all scales, enhanced by the combined effects 
of momentum conservation and non-uniformity. Also, 
the separation between fluctuations and mean flow is 
quite nontrivial. Because if the mean flow is an aver- 
age of a collection of particles, the correlations between 
the fluctuations of velocities can be big if the collection 
is big. 

An important issue is the existence of a universal de- 
scription, which is not common for nonequilibrium sys- 
tems. The separation of the dynamics into motion in the 
center of mass frame and the motion of the center of mass 
itself is quite suggestive. 

The motion of the center of mass can not be univer- 
sal. Momentum conservation decides that the velocity of 
the center of mass can be changed only by the interaction 
between particles and external effects. So it depends sen- 
sitively on the details of boundary conditions, as shown 
in the paper, and can not be universal. 

This is true for both elastic systems and inelastic ones. 
However, in elastic systems, every mode has the same 
strength due to equal-partition law of the energy. The 
motion of the center of mass is just one mode out of Nd 
modes, and its effect is negligible for a macroscopic sys- 
tem. In a dissipative system, on the other hand, being 
the only conserved mode, it can dominate over all other 
modes. Consequently, a universal description does not 
exist for the dynamics as a whole. 

Still, if we look at the other iV — 1 modes which are per- 
pendicular to this non-universal mode, we may discover 
some universal features. The independence of the motion 
in the center of mass frame on the boundary conditions 
and system sizes is a hint that this part of the dynamics 
may be universal. Further study is being carried out. 

The thin pipe model used here simplifies greatly both 
the numerical and analytical calculations. The low den- 
sity version of it may not have higher dimensional analo- 
gies, where the sequence of particles is necessarily broken. 
However, the high density version can be modified for a 
higher dimension situation, where the sequence can be 
kept. 



V. CONCLUSION 
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